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A new efficient numerical algorithm for interacting fermion systems is proposed and examined 
in detail. The ground state is expressed approximately by a linear combination of numerically 
chosen basis states in a truncated Hilbert space. Two procedures lead to a better approximation. 
The first is a numerical renormalization, which optimizes the chosen basis and projects onto the 
ground state within the fixed dimension, L, of the Hilbert space. The second is an increase of 
the dimension of the truncated Hilbert space, which enables the linear combination to converge 
to a better approximation. The extrapolation L — > oo after the convergence removes the 
approximation error systematically. This algorithm does not suffer from the negative sign 
problem and can be applied to systems in any spatial dimension and arbitrary lattice structure. 
The efficiency is tested and the implementation explained for two-dimensional Hubbard models 
where Slater determinants are employed as chosen basis. Our results with less than 400 chosen 
basis indicate good accuracy within the errorbar of the best available results as those of the 
quantum Monte Carlo for energy and other physical quantities. 

KEYWORDS: quantum simulation, strongly correlated electron systems, Hubbard model, numerical renormal- 
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§1. Introduction 

In these decades, many numerical algorithms, such as exact diagonalizations, quantum Monte 
Carlo and density matrix renormalization group(DMRG), were proposed and were applied to many 
strongly correlated electron systems. Though the exact diagonalization method is the most straight- 
forward one, the system size it can treat is smaller than that other methods can treat. Quantum 
Monte Carlo(QMC) method is a powerful technique for correlated electron systems and has been 
applied to various systems.^' In some systems such as fermion systems and frustrated spin sys- 
tems, however, the QMC is known to suffer from the negative sign problem. Namely, the cancel- 
lation of positive and negative Monte Carlo samples occurs and makes it practically impossible 
to estimate physical values in the presence of statistical and round-off errors. DMRG^) IS a very 
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powerful numerical renormalization method which does not suffer from any sign problem. However, 
because of the spatial renormalization process, DMRG is known to be applied efficiently only to 
one-dimensional configurations. 

Path-integral renormalization group algorithm(PIRG) has been proposed^^ as a new numerical 
algorithm for studying the ground state properties of strongly correlated fermion systems. The 
crucial point is that PIRG does not suffer from the negative sign problem and can be applied to 
any type of Hamiltonian in any dimension. The approximate ground state wavefunction is filtered 
out after numerical renormalization process in the path-integral formalism and expressed from the 
optimized linear combination of basis states in truncated Hilbert space. Since the explicit form of 
the wavefunction is obtained, the variational principle is satisfied and the method does not lead 
to the sign problem in contrast to QMC. Compared to DMRG, the numerical renormalization is 
done to the imaginary time direction irrespective of the spatial dimensionality of the system. Since 
the wavefunction is numerically given by using simple basis representation, it makes it easier to 
compute physical properties. 

In Chapter 2, we will explain the whole ideas and procedures of PIRG, such as renormalization 
group method, truncated Hilbert subspace and extrapolation procedure, from the viewpoint of 
physical implications and compare PIRG with QMC and DMRG. 

In Chapter 3, we will discuss the PIRG procedure from the viewpoint of implementation. There 
are some ideas and devices to make the PIRG calculation faster. We discuss the computation time 
and necessary computer memory. We review the structure of the PIRG procedure to discuss the 
parallelization of PIRG. 

In Chapter 4, we test PIRG efficiency by applying this method to the two-dimensional Hubbard 
model with nearest-neighbor transfers and compare the results with those of exact diagonalization 
and QMC. We show that the PIRG gives the results with accuracy of around three digits for energy 
and around two digits for the momentum distributions and the spin correlations. 

§2. Path-Integral Renormalization Group Algorithm 

2.1 The whole procedures of PIRG 
PIRG method consists of the following procedures. 

• Initial state. 

— It is possible to use any kind of initial state. Generally, closer to the ground state, better 
as an initial state and it is possible to use a linear combination of chosen basis states 
^initial) = Ef=i Ci|0i) or to use a single basis state 

• Projection. 

— The projection operator is introduced by exp {—tH) or by 1 — tH to achieve the ground 
state. This projection expands the stored Hilbert subspace and thus the number of basis 
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states increases. 

• Truncation. 

- Generally, the projection makes the number of basis states increase larger than our computer 
memory, or even if they can be stored, it is impossible to deal with them in the limited 
computer time. It is necessary to select important states and to keep the dimension of the 
stored Hilbert subspace L. 

• Iteration. 

— By the projection process using proper projection operator and the truncation process, the 
numerical renormalization is achieved and the stored L dimensional subspace approaches the 
ground state. It is necessary to repeat this renormalization process until the lowest energy 
of the truncated Hilbert space converges to the lowest energy under the condition that the 
dimension of subspace is L. 

Empirically, to achieve the best approximate ground state under the allowed dimension L, it is 
necessary to start PIRG using a good L = 1 state. In fact, the best result at L = 1 represents 
nothing but the optimized Hartree-Fock result. We start PIRG from L = 1 and make L larger step 
by step, using the previous converged state as the initial state at next L-dimensional subspace. 
In addition to that, it is important for a faster convergence to iterate the renormalization process 
sufficiently at small L. 

2.2 Renormalization group methods 
2.2.1 Renormalization direction 

The basic idea of the renormalization group methods is to keep the relevant information of a 
system and integrate out the irrelevant one. PIRG can be compared to the infinite-size system 
DMRG.^^ In the infinite-size system DMRG, two single sites are inserted between the system-part 
and the environment-part to enlarge these two parts as shown in Fig.l. In this process, using the 
exactly represented single-site Hamiltonians, the stored states information increases temporarily 
before the truncation process and the relevant information of the ground state is stored. This 
process is allowed only for one-dimensional system because the dimension of the Hamiltonian matrix 
of two enlarged part, or the dimension of the temporarily expanded Hilbert subspace, is too large 
to treat for boundaries of multi-dimensional systems. 

In contrast to DMRG, the numerical renormalization is performed in the imaginary time direction 
in PIRG. More concretely, the ground state \il)g) is obtained by 

\^9)=}]^^M-rH\\<t>o), (2.1) 

where \^q) is an initial state. Although the finite temprature DMRG^~^) has a similarity to this 
renormalization in the imaginary time direction, it can treat only the one-dimensional systems again 
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because the transfer matrix has to be utihzed. Generally, the state {ipg) can be represented as a 
linear combination of arbitrary basis states. However usually, the projection lim^_»oo exp[— riJ] can 
not be performed in one operation. Following the Feynman's path-integral formalism, by taking 
sufficiently small At, the projection procedure may, for example, be given as, 

IV'n) = lim lim (exp[-ATi?])"|(^o)- (2.2) 

^ Ar— »0"— ♦oo 

In quantum Monte Carlo methods, Eq.(2.2) is usually implemented at a fixed large n and sufficiently 
small Ar which are enough to filter out the ground state from the initial state \(f)o)- In the PIRG 
method, as in DMRG approach, two identical new projection operators exp[— Ar^f] are inserted 
between the stored states |V's). Namely, one operates to the right state lips) while the other operates 
to the left state (V's|- When the stored state |^s) is obtained after n projection steps from the initial 
state |0o)) the expectation value of a physical quantity (A) is represented by the initial state as 

{iPs\A\i^s) ^ (./>o|exp[-ATg]"iexp[-Arij-]"|V>o) 
(V'slV's) ((/>o|exp[-AT^]"exp[-Arfl']"|V'o) 
Here, '■<=' indicates the truncation procedure in PIRG. After the step of the PIRG procedure, in 
other word, after the projection operators exp[— Ar-ff] are inserted, this expression is changed into 
the following form. 

('(/'s I cxp [- Ar^] A exp [- Ar^] I V-s ) 
(^,| exp[-AT^] e^p[-ATH]\i)s) 
(0o|exp[-ArF]"+^ A exp[-Ari7]"+i |V'o) 
((/)o|exp[-Ari?]"+i exp[-Ari?]"+i|V'o) 

Because the dimension of the Hilbert subspace increases in this projection process, the truncation 
process is necessary similarly to DMRG method. By this projection and truncation process, PIRG 
method achieves the numerical renormalization in the imaginary time direction which is represented 
in the path integral formalism. 

2.2.2 Path integral formalism, for the Hubbard model in the Slater determinant basis 

The formalism given below shares the similarity to that of the auxiliary field Monte Carlo 
method. In this paper, we use the following Hubbard model Hamiltonian: 

H = Hk + Hjj, 

Hk = - J2 + h.c)j , 

<i,j>,a 

Hu = Uj2 (^iT - ^) {nil - ^) 

= Uj2rii^riii-u(M-^Y (2.3) 
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Fig. 1. Graphical representation of the one step of the infinite-size system DMRG. ID-system is divided into two 
parts, the system and the environment. Each of them are enlarged by one site and in this process the Hilbert space 
is also enlarged by the exactly represented one-site Hamiltonian. 



where i and j represent the lattice points, cj^ (cj^) the creation (annihilation) operator of an electron 
with spin a on the i-th site, n^o- = cl^Cia, tij the transfer integral between the i-th site and the j-th 
site, U the on-site Coulomb interaction, M = J2i,a ''^io- ^^'^ ^ the number of the lattice sites. 
The projecting operator can be divided into the kinetic and the interaction terms approximately. 
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In this paper, we assume that the basis states are Slater determinants. Hereafter, wc use the 
notation |0o-) to represent a Slater determinant with spin a and |0) = (g) \(f>i)- Since a Slater 
determinant is a single particle state, the projection of single-body operator only changes a Slater 
determinant to other single Slater determinant. Namely, 



exp 



-AHk 



(2.5) 



Though the interaction term is a many-body term, the projection term of it can be transformed 
into the sum of two single operator projections by using the Stratonovich variable s: 



exp[-ATUnm'[nmi] 
- exp [a{s)nm^] exp [Q;(-s)n„J , 

^ s=±l 



where 



a{s) = 2as 



AtU 



, _i , /AtU\ 
tanh W tanh ( — - — I . 



(2.6) 

(2.7) 
(2.8) 



As a consequence. 



exp [-ATUnmirimi] 



(ir+) + ir-)) 



(2.9) 
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exp 



(2.10) 



i=l 



In this way, the projection of the local interaction term, exp [— Art/n^-i-n^j^], changes a Slater 
determinant to the sum of two Slater determinants. After the projection of the AT-sites interaction 



term exp 



-At Hi 



u 



, an original single Slater determinant expands to the sum over 2^ Slater deter- 
minants. Though we use Slater determinants and a projection exp[— Ar^f], a similar dimensional 
expansion occurs even if we use other kind of basis states such as the site representation and a 
projection 1 — H. This is the reason why QMC sampling or PIRG truncation is necessary. 

2.3 Truncated Hilbert space 
2.3.1 Variational principle 

The expectation value {)g of a physical variable A in the ground state I'^p) can be calculated 
with an arbitrary complete set of basis states as follows: 

^ Hub ert 



(A 



E 



(2.11) 



However, the sum in the above equation is practically impossible in general for easily available 
basis states because the number of them, D Hubert, is usually comparable to the dimension of the 
whole Hilbert space. QMC methods deal with this problem by sampling the numerator and the 
denominator separately: 

/ Probability Sign Sample 



(a,6) 



E 



X 



(0a|^6) 
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I 



E 



(cd) 



\ 



(2.12) 



E.^ 



Dh 



V 
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Sign 



where || || represents the absolute value and Ng, the number of Monte Carlo samples. Usually, these 
sampling processes arc taken for an identical set (c, d) = (a, b) in the Metropolis algorithm. 

Because the diagonal elements {(l)a\<Pa) are not necessarily contained in Monte Carlo samples and 
{4'a\(t>b) in th<3 sum are taken only partially by sampling, QMC does not satisfy the variational 
principle in the strict sense, although the deviation should be within the range of statistical error. 
Several thousands samples are practically taken, to reach the accTiracy with a relative error less 
than a few pcrccnts. In case of the system with the sign problem, however, it is difficult to achieve 
the same accuracy by this sampling process. 
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On the contrary, the variational principle is satisfied and the sign problem is absent in PIRG 
because the ground state is represented approximately as a linear combination of arbitrary basis 
states. 

L 

\i^g)^H)=Y.'^Mi)- (2-13) 

i=l 

In Eq.(2.13), are optimal basis states in the whole Hilbert space. For example, orthogonal basis 
states such as site-represented ones or non-orthogonal basis states such as Slater determinants can 
be used. In any case, to take the sum in Eq.(2.11) within the allowed computation time, these 
basis states should be simple ones and the number of them, L, be small. Of course, the relation 
between L and the truncation error may depend on the choice of basis states. Although we discuss 
the relation next, we do not know it clearly now. Empirically, as we show the results in Chapter 4, 
with the choice of the Slater determinant bases, limited and tractable basis states, such as hundreds 
basis states, can reach the ground state with the relative error less than a few percents for most of 
physical quantities of our interest even when the Hilbert space is enormously large. 

From theoretical and practical viewpoints, the relation between the error of estimated properties 
and the truncated Hilbert space dimensions L is very important. It is known that the relative 
error decreases exponentially in DMRG^^ when the number of states kept, or the dimension of the 
subspace, increases. This makes DMRG powerful. In DMRG, because the states are not treated 
explicitly, there is no restriction on the choice of bases. Namely, although both PIRG and DMRG 
are the methods to make the truncated Hilbert space converge to the ground state, the restrictions 
posed on the choice of basis states are different in two methods. 
2.3.2 The lowest energy in the truncated Hilbert space 

As shown in Eq.(2.13), the stored approximate ground state is represented as, 

L 

1^) = ^Wi\(t)i) 
1=1 

Because the basis states which constitute the stored Hilbert subspace are chosen so as to give 
the lowest energy within the allowed dimensions of the Hilbert subspace, it is necessary to optimize 
the coefficient Wi in every truncation process. The energy E of the state jV') follows the equation, 

where, 



E = t:^/=^;;:'^ ' ^ (2.14) 



[F\ij = (2.15) 
The lowest energy E in this subspace is the lowest eigenvalue of this subspace. Then, to obtain 
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the lowest energy and its state, we solve the following generalized eigenvalue problem. 

L L 

Y,[H]ijWi = EY,[F]iiWj. (2.16) 

i=i i=i 

2.4 Extrapolation on the dimension of the stored subspace to the dimension of the whole Hilhert 

space 

2.4-1 Extrapolation of energy 

We need the extrapolation procedure to large L to estimate the systematic deviation, because 
the exact value is achieved only when L becomes the dimension of the whole Hilbert space. In 
practice it is difficult to analyze the relation between the dimension L and the systematic deviation 
of expectation values. In addition to that, the relation may depend on the choice of the bases. 
Therefore, we do not use the argument L in the extrapolation function and introduce another 
extrapolation procedure. The important point is that the extrapolation procedure gives more 
accurate estimate. Wc, however, note that after this extrapolation, the variational principle may 
not necessarily be satisfied because the obtained energy could be lower than the exact one within the 
extrapolation error. We define the difi^erence between the ground state energy and the expectation 
value in a given subspace as 

5E = {H) - {H)g (2.17) 

It can be shown^) that the difference SE vanishes linearly as a function of the energy variance AE 
defined by 

AE = ^^'^ 7 (2.18) 

We summarize the proof of this in the following way: First we put the approximate ground state 
IV') as 

\i;)=c\i;g) + d\i;e) (2.19) 

c2+d2 = l 



where \ipg) and |V'e) are orthonormalized states. We also define 



{H)e - {H)g 



Do = 



{H? 



a 



where 



'9 

D3 = - ' \ ^ , (2.20) 

m 



^ = M^IM, (i). = (2.21) 
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In this notation, 



{H) = {H)g+S{H)yDl 

{H^) = {H)l + <f{H)lD^ 
{H') = {H)l + d'{H)lD, 



SE = cPEqDi 

AE = (f {D2 - 2Di) + d'^Di (3£)i - 2D2) + O {S) 

(2.22) 

where {A) = {ip\A\'il)) / {'tp\'ip) and Eq = {H)g. When the stored state is a good approximation 
of the ground state, the coefficient d is expected to be small. Then up to O (d^) , 

5E oc AE (2.23) 

is satisfied. This is the simplest extrapolation procedure wc introduce. 

We may introduce other scries of extrapolation procedure. They arc more time consuming but 
can be more accurate. As wc discussed in the rcnormalization procedure, there arc operators that 
can make a state closer to the ground state. Here wc take this kind of projection operator R. In this 
case, for more accurate extrapolation, 6E and AE in the above simplest extrapolation procedure 
should be calculated on the state R\ip). Namely, in the above equations for any operators A, 

(i> = (2.24) 

may give a better estimate than that from (iplAlip) / {iplip) . In the rcnormalization procedure, we 
have to consider the restriction of the Hilbcrt subspace dimension due to computer power. As 
shown in the above, however, if it is possible to calculate in the form Eq.(2.24) with the projection 
R, it is equivalent to calculate the physical value with the state R\tp), which is in a larger Hilbert 
subspace than that the state lip) belongs. 

As an example of the projection operator R, we take a function of Hamiltonian H. When 
the absolute value of the ground state energy is the largest in all the eigenvalues, operating the 
Hamiltonian H makes a state closer to the ground state. From a viewpoint of computation time, the 
number of the terms contained in the Hamiltonian of the system without any long-range interaction 
is the order of N. If we take the Hamiltonian H as the projection operator R, it is necessary to 
calculate {H'^) in AE but it is not practical because the computation time increases to the order 
of N'^. Therefore in this paper we take \/h as the projection operator. To obtain the expectation 
value of A on the state i2|(/!>), we use Eq.(2.24). In this case, Eq.(2.22) is modified to the following. 

SE ^ SEsgrt = 7-7 - (^)9- (2.25) 



9 



A£; ^ AEsgrt = ' '\;' ' ■ (2.26) 



{H^){H)-{H^f 

The difference dEsqrt is always smaller than the difference SE and the extrapolation procedure can 
be done closer to the ground state if we use SEgqrt and AEgqrt- One might consider a problem 
whether the operator exists or not. In the representation of SEgqrt and AEgqrt, however, 
does not appear but only H does and the argument for SEggrt and AEgqrt remains correct. Even 
if we ignore the projection operator R = in Eq.(2.24) and only use Eq.(2.26), the following 
equations are obtained. 

' ?>Egqrt = d^EQ{D2-D^) 

+d'^EoDi (2Di -D2) + (cZ6) 
AEgqrt = d"^ {D3 + Di- 2D2) 

{D3D1 - 2D2D3 - 2D1D2 + 3D|) + O (d^) 

(2.27) 

By the same reason as that of the simple extrapolation procedure, the above equation leads to 

SEgqrt OC AEgqrt- (2.28) 

Because both extrapolations are on the energy difference and the energy variance, they can be 
plotted on the same parameter plane. As shown in the above, however, two extrapolation plots are 
on different lines. 

SE ^ AE (2.29) 

D2-2D ^ ' 

where the proportionality coefficients are different between Eq.(2.29) and Eq.(2.30). In this paper 
we use the SEgqrt — AEgqrt extrapolation to examine the accuracy of the ground state energy. 
2.4-2 Extrapolation of other physical values 

Because no commutation relation is expected between most operators of the physical quantity A 
and Hamiltonian H., 



{^\VHVH\ct>) {m\(b) 
Then we can not use the above second extrapolation procedure. From the same reason as the 

computation time for obtaining the energy, higher order projection is not practically possible within 

our accessible computer power. For the moment, we do not know the suitable projection operator R. 

Then, we use the simple extrapolation for other physical values than the energy. Some modifications 

are, however, necessary in the above discussion. The expectation value of A is 

(A) = {A)g + d^{A)gDA + 2cd{ilJe\A\i^g), (2.31) 



Da = 
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Different from the discussion on the extrapolation of the energy, in which the energy difference and 
the energy variance is proportional to ^ , the third term in (A) is proportional to d and therefore, 
the extrapolation of the expectation value of the operator A should be done on the square root 
of the energy variance. Empirically, however, the d term is so small that it is easy to do the 
extrapolation on (A) linearly to the energy variance AE, and we do so in this paper. Theoretically, 
it has been proven^) that if the operator A is a short-ranged correlation function, {ipelM'^g) 
Eq.(2.31) becomes zero in the thermodynamic limit as the following: Because \ipg) and IV'e) are 
orthonormalized states. 



(i (2.32) 
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where the latter inequality is led by Schwartz inequality and ( )g = {ipgl \ipg)- When yl is a 
short-ranged correlation function, the final term is proportional to 1/N and becomes zero in the 
thermodynamic limit. 

These extrapolation procedures on energy and other physical values are satisfied in a general 
ground state. They arc not necessarily restricted to the purpose of the application in PIRG method 
and not necessarily for studying the ground state. It is possible to generalize these extrapolation 
procedures for obtaining a more accurate eigenstate properties from a series of approximations of 
eigcnstatc. 

2.4-3 Comparison with previous methods 

Truncations of the Hilbert space were used in many numerical method. Although DMRG is 
one of the most remarkable approaches, many other algorithms applicable to higher-dimensional 
systems have also been proposed. 

Exact diagonalization method combined with the truncation of the Hilbert space^"^^^ can treat 
lagrcr systems than a normal exact diagonalization method can do. The efficiency of the usage 
of basis state, however, is not sufficient partly because the basis states are in site-representation. 
Therefore the convergence of the energy as a function of the dimension of Hibert subspace is much 
worse than that of PIRG which uses Slater determinants for basis states. 

On the other hand, the quantum Monte Carlo method is also combined with the truncation 
of the Hilbert space and applied to fermion systems"*^^) and boson systems. ""^^^ Although in these 
methods, candidate basis states arc generated by Monte Carlo process of the whole Hamiltonian, 
in our experience, it is crucial for the lowering of the energy to make the acceptance of them higher 
by sequentially generating local, or one-site, projection. Renormalization and projection through 
the interaction term is achieved more cfficeintly by the local algorithmas we discuss in § 3.2.1. 
Several devices to get out of local minima and to realize global minimum are important in models 
of condensed matter systems. 
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A large difference between tlie above numerical algorithm and PIRG is the presece of the combined 
extrapolation procedure. It can make errors smaller systematically. 

§3. On the implementation of PIRG 

3.1 Matrix representation 

The basic ideas of the following discussion are published. At first, we explain the expression of 
states and expectation values. A Slater determinant \(j)a) is represented as an A'" x M matrix [(t)c\- 

M / N 



M / N \ 

'a) = li \12i^-h4a] \0) (3.1) 
7=1 \i=l / 



The inner product of two Slater determinants is 

{4>aa\4>ba) = det (* [(j>aa] [(l>ba]) ■ (3.2) 

(Sec Appendix). The matrix elements for the Hamiltonian and other operators are calculated from 
single-particle Green's function using the Wick's theorem, because a Slater determinant is a single- 
particle state. A single-particle Green's function [Gf^], which is represented as an A^ x A^ matrix, is 
calculated by 



ab^ . . ^ {(l)aa\4a(^ja\4>ba) 
'■^ {(l>aa\(pba) 



M M 

j:j:[<f>baU\9f]j<i>aa]jl (3.3) 



k=l 1=1 

where 

[9f] = {'[<l>aama]y'. (3.4) 

The procedure to use the Wick's theorem is explained in the literature.^) Since the stored state in 
PIRG is a linear combination of Slater determinants, the expectation value of physical variables 
can be calculated from one-particle Green's functions between all the stored Slater determinants 
by using the Wick's theorem. 

Next we explain the projection processes following the above representation. The kinetic term 
projection, Eq.(2.5), is expressed by using an N x N matrices [Mq] and [K] in the following way: 

N 

Wahj = T.i^o]ik[<t>a]kj for a =T, i (3.5) 
k=l 

where 

[Mo] = exp[K] 

[K]ij = -Artij. (3.6) 
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Prom Eq.(2.6) the local-interaction term projection, Eq.(2.9), is expressed as the following: 



/ N \ / ^ 

\k=l 
/ N 

E [^1 (-1)]^;^ 



\k=l 



N 



where Mi is an A'^ x A'' matrix, 



[Ml (s)] 



exp 

1 





2as 



J kj 



AtU 



E [^1 (1)]^;^ 



(3.7) 



for i = j, 
otherwise 



\A;=1 

for i 

i ^ m 



J = m 



(3.8) 



3.2 Computation time 

3.2.1 Whole procedures of the renormalization group method 

Empirically we find that it is better to perform the projection and truncation at the projection 
of every local-interaction term exp [— ArJ7nmt?^mJ- Hence we employ this local algorithm. Using 



Ea=l Ca\(i)a) 



this local algorithm, we summarize below how one PIRG iteration step exp 
proceeds. Note that in this paper represents \(pa]) ® l^ai); namely a direct product of Slater 
determinants for each spin. 

• Choose one basis state. 
— Choose \<t)a) from L stored basis states {|0i), 1^2); • • • , I^^'l)}, which will be operated by 



exp 



-Ari7 



• Projection by exp 



-AtH. 



com,putation time oc [LN^ + L^) 



— Perform 



exp 



-AtHu 



Hamiltonian elements between 



and 



\(pa) following Eq.(2.5) and calculate the inner products and 

\4>'a) 



{\(f>l), Ih), \(l>a-l), \(f>a+l), IH)} ■ 

— > computation time oc LN^ 

Truncation is performed by comparing the lowest energies obtained from Eq.(2.16) in two 
subspaces which consist of the following two sets of L basis states: 

{\4'l), \4'2), \4'a^l), 10a), \4'a+l), \4>l)} 

{\4'l), \4'2), {(Pa-l), \(l)'a)i \4'a+l), 

and by employing one of these two basis states set which gives the lower energy. In other 
words, take \(pa) or {(p'^) to be a next basis state |0°). 
— > computation time oc 
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• Projection by exp —AtHu . 

=^ computation time oc (LN^ + L^N) 

— Perform 

I m-") + IC-)) = exp [-ArUnm^rimi] IC^ 
following Eq.(2.9) and calculate the inner products and Hamiltonian elements between 

or IC") 

and 

{|0l), 102), ■ • • , \(f>a-l), \(l>a+l), ■ ■ IH)} ■ 

computaiton time oc LN'^ 

— Truncation is performed by comparing the lowest energies obtained from Eq.(2.16) in three 
subspaces which consist of the following three sets of L basis states: 

\^2), . . . , \cpa-l), \<t>T-'), \<Pa+l), 
{|</.l),|(/>2),...,|0a-l),|0™+),|</'a+l),...,|'/'L)} 
{|</.l), |<^2), . . . , |0a-l), \^a+l), ■ ■ ■ , 

and employing one of these three basis states set which gives the lowest energy. The chosen 
basis state \(t)'^~^), 10^"*") or|0^~) is taken to be a next basis state |(?!>^). 
— > computation time oc 

— Repeat the above local interaction projection for all the sites m = 1,2,...,N. 

• Repeat the above exp —AtH \(f)a) projection for all basis states a = 1, 2, . . . , L. 

=^ The total computation time oc L^N^ + L'^N (3-9) 

Though the basic computation time of PIRG is listed above, an efficient convergence procedure 
is important for reducing the computation time. There may appear many states with local minima 
structure in energies in the PIRG convergence process. Even if the convergence is not perfect, 
which is actually the case in most of our experience, the extrapolation procedure can be performed 
for the ground state properties. However, the worse converged state gives the value with larger 
error. Because the extrapolation itself using the energy variance is the formalism to obtain values 
of an cigcnstatc of the Hamiltonian, in case of worse convergence than a limit, the extrapolation 
procedure could give the value not for the ground state but for an excited state. Therefore it is 
important to improve PIRG convergence process by a combination with some existing methods to 
avoid occurrence of trapping in local minima. At present, we have empirically learned that the 
sufficient convergence at small subspace is crucial at the early stage of the process of extending L. In 
this paper, we have numerically realized convergence to the state of the Hartree-Fock approximation 
at L = 1 and iterated hundreds projections at small L, for example for L less than 10. However 
more systematic methods are desired and are left for future studies. 
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3.2.2 Details on reduction of computation time 

Since inner products and Hamiltonian elements are calculated from determinants and inverse 
matrices, the computation time for them is usually proportional to N^. In the process of local- 
interaction projection, however, it can be reduced to N"^. 

Here, to explain the procedure to reduce the computation time from to N^, we simply consider 
the state J2k=i [-^i (^)]jfe [<?^ao-]fej Eq.(3.7) for the following discussion. When the basis state of the 
a-th Slater determinant with spin a is updated by the projection of the m-th site local interaction 
as. 



[<Pa 



N 



\ik Vraalkj 



(3.10) 



the inner products change as, 

{(l>aa\(l>ho) -■ 

and the Green's functions change as, 

(<5 (s) Sim + 1) 



{l + 5{s) 



ab 



ab 



) for 6 7^ a 



Wl^feff) for6 = o 



(3.11) 



g: 




'Gf 


. H 

im 


s) 


Gf 


mj 


1 + 


Gf 


5{s) 

mm 



Gcr 


. H 

im 


s) 


Gcr 


mj 


1 + 




6{s) 

mm 



for 6 7^ a 



for b = a 



(3.12) 



where 



G„ 



{S (s) Sim + 1) 



g: 



ab 



G' 



ah 



6{s) 



G 



ab 



mj 



1 + 



g: 



ab 



d{s) 



S (s) = exp 



2as 



AtU' 



1. 



(3.13) 
(3.14) 



Here i,j is from 1 to and 6im is the Kronecker's 6 symbol. In this way the computation time 
for local interaction projection is reduced to that proportional to A'^^. The similar reduction in the 
computation time is explained in detail in the literature."^) 
3.2.3 Devices on extrapolation procedure 

For the extrapolation process after the PIRG convergence, it is necessary to calculate (H), {n), 
{H'^) and expectation values of other operators. Because the Hamiltonian does not have any 
long-range interaction term, the computation time for (if") is proportional to A^" except for the 
computation time for single particle Green's function N'^. The coefficient of N'^ is, however, large 
and it is useful to decrease this computation time by the following procedure. 

Because a converged state by PIRG is a linear combination of Slater determinants and the kinetic 
term projection exp[— AriJ/j] does not increase the number of Slater determinants, it is easier to 
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calculate the right hand side of the following equation than to calculate the left hand side from 
single particle Green's functions in the following equation: 



AHM) 1 / (MAm (V'l^exp 



(VIV') At \ {^^P) (V'lV) 



+ 0(At), (3.15) 



where we take sufficiently small At. In this way, the computation time of the term, such as H^, 
can be reduced. 

3.3 Required memory 

The largest share of the memory is exhausted for the storage of the elements in all the elements 
of the Slater determinants in the basis states and also for that of the Green's function. Here we 
assume scalar information takes 8byte and the number of basis states L is 500 to roughly estimate 
the required maximum memory size. Slater determinant is an x M matrix for each spin, and 
then near half filling of the Hubbard model, this is comparable to N"^ scalar elements. As we see in 
(3.3), Green's function is represented by an x A'' matrix for each spin, and then this contains 2N'^ 
scalars. If we store all the Green's function data among the L Slater determinants, about 83Gbyte 
is necessary for 12 x 12 system and it cannot be easily stored in our available computer. For this 
reason, we store only the Green's function between the Slater determinant just on the operation of 
exp[—ATH] and the others. Then the necessary memory for the state and Green's function is the 
order of 3LA'^^ x 8byte, which is about 250Mbyte for 12 x 12 system. 

3.4 Parallelization 

Parallelization of a code is a promising way to improve the performance of the computation by 
dividing a large set of calculations into several smaller pieces and execute them independently of 
each other. We have tried to parallelize the code on the distributed memory system using MPI. As 
shown in Eq.(3.9), the computation time for a single slice projection exp[— AtJ?] is proportional 
to L x (LN^ + L^N) . Because the first factor L refers to the iteration on all the basis states which 
are related to each other in the evaluation of the energy, it is difficult to parallelize this process. 
Then we consider the parallelization of the part LN^ + L^N. 

The first term LN^ is related to the calculation of inner products and Green's functions between 
a basis state on the operation and the other basis states. Each calculation is independent and it can 
be parallelized. For this parallelization, the memory for each state and each Green's function have 
to be distributed over each processor. For example, the 1-st state and the 1-st Green's function are 
on the 1-st processor memory, and the 2-nd state and the 2-nd Green's function are on the 2-nd 
processor memory, etc.. 

The second term L^N is related to the iteration of the calculation of the lowest energy of the 
stored Hilbert subspace by Eq.(2.16) for each local projection. The matrix related to this generalized 
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eigenvalue problem is L x L. In practice, L is the order of hundreds and it is difficult to parallelize 
efficiently the eigenvalue problem of such a small matrix. However by performing the projection 
exp[— Arff] for some different choices of At in parallel and employing the result which gives the 
lowest energy among the choices of At, the convergence becomes faster and computation time for 
Eq.(2.16) does not increase. In this way, we can reduce the total computation time for convergence 
by parallelization. 

§4. Evaluation of PIRG on Hubbard model 

4.1 Model 

We apply PIRG to the Hubbard model Hamiltonian (2.3) on a two-dimensional square lattice 
with nearest-neighbor transfers t. 

{t = 1.0 ii are the nearest neighbor sites 
otherwise 

We take t = 1.0 as the energy scale. Since these models have particle-hole symmetry at half filling, 
the quantum Monte Carlo does not suffer from the sign problem at half filling. By comparing 
with the QMC and exact diagonalization results in the literature, we show PIRG results on the 
energy, the momentum distribution and the equal-time spin correlations and discuss the accuracy 
and efficiency of PIRG on these models. 

4-2 Comparison of energy 

4.2.1 PIRG result 

First, we show the converged results by PIRG before the extrapolation procedure and discuss 
the relative error. In §2 we refer to the DMRG method in which the relative error decreases 
exponentially as a function of the dimension of the stored Hilbert subspace. Because there is a 
restriction on the choice of basis in PIRG, the relative error dependence on the dimension of the 
stored Hilbert subspace is different between DMRG and PIRG. Figure 2 shows the relative error 
and Figs. 3 and 4 show the relative difference between PIRG and QMC results. In case of a small 
system such as 6 x 2 Hubbard model, the relative error is smaller than 1 percent. For larger systems 
near half filling, the relative error is larger but less than a few percent independent of the system 
size. This feature holds both at half filling and near half filUng. Note that the results at L = 1 
correspond to the Hartree-Fock estimates. 

4.2.2 Extrapolation results 

The above results are improved to more accurate estimates by the extrapolation procedure. Here 
we show the simple extrapolation results E with the state |^) and the extrapolation results Egqrt 
with the state ^/h\4>) for the same 6x6 Hubbard model by taking L up to 256. Both extrapolation 
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Fig. 2. Relative error SE/\E\ in the ground state energy for the 6x2 Hubbard model with 5 up 5 down electrons 
and the fully periodic boundary condition at U/t = 4.0. The reference ground state energy is estimated from the 
exact diagonalization. 
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Fig. 3. Relative difTerence 5E/\E\ in the ground state energy for the Hubbard models at half filling with the fully 
periodic boundary condition at U/t = 4.0. The reference ground state energy is estimated from QMC.^^ 



procedures for two models are shown in Figs. 5 and 6. Two extrapolated results should meet at 
the same value in principle, although there is a small difference, which may be thought to be PIRG 
error. Because we do not know which one is closer to the exact value and empirically it is better for 
the linear function fitting to use the value obtained from the state hereafter we employ the 

extrapolation using \/H\(p). The relative errors and differences after the extrapolation are shown 
in Table I. For most of the systems, PIRG can give results of the ground state energy with less 
than 0.3 percent relative error or difference from QMC results. 
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Fig. 4. Relative difference 5E/\E\ in the ground state energy for the Hubbard models with two hole doped from 
half filling and fully periodic boundary condition at U/t = 4.0. The reference ground state energy is estimated 
from 
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Fig. 5. Extrapolation of the energy to the zero energy variance for a 6 x 6 Hubbard model, 17 up 17 down electrons 
with the fully periodic boundary condition at U/t = 4.0 



4-3 Comparison for other physical quantities 

Next, we evaluate the equal-time spin correlations and the momentum distribution on 6 x 2 
Hubbard model, 5 up 5 down electrons with the fully periodic boundary condition atU/t = 4.0. In 
our study, the equal-time spin correlations in the momentum space is calculated from, 

^ (9) = ^ E iSiSj) eMR^-^^) (4.2) 
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Fig. 6. Extrapolation of the energy to the zero energy variance for a 6 x 6 Hubbard model, 18 up 18 down electrons 
with the fully periodic boundary condition at U/t = 4.0 



Table I. relative error and difference of the ground state energy after extrapolation. 



system 


PIRG results 


exact diagonalization and Monte Carlo results 


relative difference or error 


6x2, 5T,5i 


-25.6999 ± 0.0005 


-25.6952 


0.00018 


6x6, 17T,17i 


-65.12 ± 0.02 


-65.30 ± 0.04 


0.0028 


6 X 6, 18 T, 18 i 


-66.92 ± 0.04 


-66.96 ± 0.07 


0.00060 


8 X 8, 31 T,31 i 


-117.8 ± 0.1 


-117.70 ± 0.06 


0.00085 


8 X 8, 32 T,32 i 


-119.4 ± 0.1 


-119.23 ± 0.06 


0.0014 



where Si is the spin of the i-th site and each element of the spin is calculated from 

St = \ {st + S-) = \ (cj^ci + 4^c,t) 

St = \{n,^-n,^). (4.3) 

All the extrapolation behaviors of the equal-time spin correlations are shown in Figs. 7 and 8. The 
results after the extrapolation procedure are shown in Table II. Here we take the lattice constant 
to be unity. 
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Fig. 7. The extrapolation of the equal-time spin correlations S {kx, ky) to zero energy variance for the 6x2 Hubbard 
model with 5 up 5 down electrons with the fully periodic boundary condition dXU /t = 4.0. The wavenumber {kx , ky) 
for each symbol is given in the inset. 
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Fig. 8. The extrapolation of the equal-time spin correlations S {kx, ky) to zero energy variance for the 6x2 Hubbard 
model with 5 up 5 down electrons with the fully periodic boundary condition at U/t = 4.0. The wavenumber {kx, ky) 
for each symbol is given in the inset. 



The momentum distribution is calculated from, 

" (9) = ^ E (sV^T + e<^-R^) (4.4) 

where Ri is the vector representing the place of the z-th site. All the extrapolation behaviors of 
the momentum distribution are shown in Figs. 9 and 10. The comparison of the results after the 
extrapolation is shown in Table III. 

Among the equal-time spin correlations and the momentum distribution, some of them have 
relatively large errors of about a few percents. Most of the relative errors are, however, less than 1 
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Table II. The equal-time spin correlations S {kx, ky) for the 6x2 Hubbard model with 5 up 5 down electrons with 
the fully periodic boundary condition at U jt = 4.0 



{kx-, ky) 


PIRG 


exact diagonalization 


relative error 


(1,0) 


0.0482 


0.0488 


0.012 


(2,0) 


0.0454 


0.0458 


0.0044 


(3,0) 


0.0451 


0.0457 


0.013 


(0,1) 


0.281 


0.277 


0.014 


(1,1) 


0.283 


0.281 


0.0071 


(2,1) 


0.284 


0.286 


0.0070 


(3,1) 


0.278 


0.279 


0.0036 



percent. 
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Fig. 9. The extrapolation of the momentum distribution n{kx,ky) to zero energy variance for the 6x2 Hubbard 
model with 5 up 5 down electrons with the fully periodic boundary condition at U/t = 4.0. The wavenumber 
{kx, ky) for each symbol is given in the inset. 



§5. Summary 

Path-integral renormalization group(PIRG) is a numerical algorithm for studying the ground 
state properties. The process filtering out the ground state \ipg) is performed in the imaginary time 
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Fig. 10. The extrapolation of the momentum distribution n {kx, ky) to zero energy variance for the 6x2 Hubbard 
model with 5 up 5 down electrons with the fully periodic boundary condition at U/t = 4.0. The wavenumber 
{kx,ky) for each symbol is given in the inset. 



Table III. The momentum distribution n {kx,ky) for the 6x2 Hubbard model with 5 up 5 down electrons with the 
fully periodic boundary condition at U/t = 4.0 



{kx, ky) 


PIRG 


exact diagonalization 


relative error 


(0,0) 


0.9678 


0.9681 


0.00031 


(1,0) 


0.9604 


0.9601 


0.00031 


(2,0) 


0.9288 


0.9281 


0.00075 


(3,0) 


0.02067 


0.01850 


0.092 


(0,1) 


0.06515 


0.06806 


0.044 


(1,1) 


0.04504 


0.04513 


0.0017 


(2,1) 


0.02821 


0.02787 


0.012 


(3,1) 


0.02294 


0.02281 


0.0057 



direction as, 

Therefore its formalism can be applied to any kind of systems and there is no restriction on the 
spatial dimension of the system. In this path-integral formalism, the ground state is represented 
by chosen basis states |^), 
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By the numerical renormalization, relevant basis states are selected and irrelevant basis states are 
projected out. This makes it possible to calculate the approximate ground state directly as an 
optimized linear combination of chosen basis states: 

L 
i=l 

Therefore there is no negative sign problem even in frustrated systems. In this way, PIRG can be 
applied to the systems which can not be treated by existing algorithms such as the quantum Monte 
Carlo method or the density matrix renormalization group. 

Because the converged state by PIRG is an approximate ground state under the restriction on 
the number of the basis states, the exact ground state can be achieved by the extrapolation of 
the number of states L to the dimension of the whole Hilbert space. We have shown the general 
extrapolation procedure in this paper: 

{H) - {H)g oc AE 

where AE is the energy variance, 

( )g, the expectation value in the true ground state and ( ), the expectation value in an approx- 
imate ground state. This relation holds for sufficiently converged approximate state. We confirm 
that the results of more than a hundred Slater determinants follow the above relation and can 
be used for a linear extrapolation in the Hubbard model. On the physical quantity A, a similar 
relation holds in most cases. 

(i) - {A)g (X AE. 

For short-ranged correlation functions, the linearity holds at the same level as the energy in the 
thermodynamic limit. By these extrapolation procedure, more accurate results are obtained while 
the variational principle is not strictly satisfied after the extrapolations. 

We compare the expectation values of the energy, the momentum distribution and the equal- 
time spin correlations with those of exact diagonalization on 6 x 2 lattice systems and with those of 
quantum Monte Carlo on larger systems. For the momentum distribution and the equal-time spin 
correlations, the relative errors and differences from QMC are less than a few pcrccnts. Especially 
on the energy, the relative errors and differences from QMC have three digits accuracy. We confirm 
that these accuracy can be achieved up to 12 x 12 systems on the square lattice. 

We have also explained the dependence of the computation time as Lp'N'^ + L'^N where N is 
the system size and L, the dimension of the stored PIRG Hilbert subspace. We refer to some 
implementation advice on PIRG and the efficiency of distributed memory parallelization on PIRG 
methods. Because in general it is difficult to parallelize efficiently the operations of matrices with 
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the size such as hundreds x hundreds or the iteration process, it is necessary to change the algorithm 
to make the convergence faster by projecting some basis states in parallel. 

In this study, we have dealt with Hubbard models using Slater determinants as PIRG basis states. 
Applications of PIRG to other systems are very interesting and promising future projects. 
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Appendix: Inner product of Slater determinants 

Here we ignore the spin degrees of freedom for simplicity. A Slater determinant \(j)a) represented 
to be an A?^ X M matrix [cpa]- 

M / N \ 

K)=n E[^«]^.-4 io). (A-i) 
j=\ \i=i / 

where a is a symbol to distinguish Slater determinants. Eq.(3.2) can be obtained as the following: 

M M / N \ f ^ \ 

= (oi n n E ['^j.. ^0 E t'^^]'^' ^^-2) 
j=ii=i \i=i J \k=i J 

nCm ( M\ M\ M , . 1 

Here n and m are permutation of M symbols: 

n=\ , (A-4) 

\^ n(l),n(2),---,n(M) ) 

the same is assumed for m; Ylm' means taking the sum over all permutations of M symbols; sgn{n) 
is the signature of the permutation n: 

{+1 for an even permutation 
; (A-5) 
— 1 for an odd permuation 

S is a set of M numbers chosen from N numbers 1,2, - ■ ■ ,N and we assume the ascendent order 
on this set S; J2g means taking the sum over all combinations of M numbers. Then the inner 
product can be transformed from Eq.(A-3) as the following: 

nCm ( Ml M 



NyM m '■ Lvi , s 

{<t>a\4>b) = ^' E |E ^5n(n) n {\.<i>a]s,j W^nop- j | (A-6) 

nCm ( M\ M ^ 



n j=i 
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M! M / N N 

^ sgn{n) \\ ^ [<i)a\j [Mmij) 
j=i \i=i J 



n 



det (* [</.„] [cl>,]) . 



(A-8) 
(A-9) 
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